//ssc install estout
//ssc install rwolf2  

global replica_dir="<add path to your working directory>"
global data_dir="$replica_dir/data"
global output_dir="$replica_dir/output"

use "$data_dir/panel_simce4b8_sned_studlvl.dta", clear


*------------------------------------------------------------------------------*

 
 
global depvars "mate lect nat "

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

la var treat25 "Full vs. Partial Bonus"
la var treat10 "Partial vs. No Bonus"

la var st_ind25 "100% vs. 60%"
la var st_ind10 "60% vs. none"

 
foreach yvar of global depvars{
	foreach i in 25 10{

		rdbwselect `yvar'_post st_ind`i'  					, all  c(0) p(2) kernel(uni) masspoints(off) vce(cluster rbd)

		gen bw_mse_`yvar'_`i'=e(h_mserd)
		gen bw_cer_`yvar'_`i'=e(h_cerrd)

		rdbwselect `yvar'_post st_ind`i' if sch_municip_pre==0 	, all c(0) p(2) kernel(uni) masspoints(off) vce(cluster rbd)

		gen bw_mse_`yvar'_`i'_pv=e(h_mserd)
		gen bw_cer_`yvar'_`i'_pv=e(h_cerrd)

		rdbwselect `yvar'_post st_ind`i' if sch_municip_pre==1 	, all c(0) p(2) kernel(uni) masspoints(off) vce(cluster rbd) 

		gen bw_mse_`yvar'_`i'_mv=e(h_mserd)
		gen bw_cer_`yvar'_`i'_mv=e(h_cerrd)

}
}

foreach i in 25 10{
	foreach y of global depvars{

	replace bw_cer_`y'_`i'=round(bw_cer_`y'_`i',0.01)
	replace bw_mse_`y'_`i'=round(bw_mse_`y'_`i',0.01)

}
}



save "$data_dir/temp_b8_bw_nocovs.dta", replace


*------------------------------------------------------------------------------*




/*
	Regression RD
*/

use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat " // soc_post 

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"
	
foreach i in 25 10{

local label_`i': variable label st_ind`i'

foreach y of global depvars{
	
		local label_`y': variable label `y'_post

		matrix define R = J(6, 9, .) // ( #of bw , #of statistics, missing )


		sum bw_cer_`y'_`i'
		local cer=round(r(min),0.1) // r(min) // 
		di `cer'

		sum bw_mse_`y'_`i'
		local mse=round(r(min),0.1) // r(min) // round(r(min),0.1)
		di `mse'
		
		local cer2=`cer'*2
		local mse2=`mse'*2
		
		local r = 1
		foreach bw in `cer'  `mse' `cer2'  `mse2'   1.5  3 {
				
		
		reg `y'_post treat`i' st_ind`i' st_ind`i'_sq inter`i' inter`i'_sq  ${indepvars} `y'_pre  if abs(st_ind`i')<= `bw', vce(cluster rbd)

			matrix list r(table)

			matrix R[`r', 1] =`bw' //round(`bw', 0.1) // name of bandwidth 
			matrix R[`r', 4] = r(table)[1,1] // rd estimate ****** used
			matrix R[`r', 5] = r(table)[2,1] //se_ of rd estimate ****** used
			matrix R[`r', 6] = r(table)[4,1] //2 * normal(-abs(R[`r', 4] / R[`r', 5])) // p-value
			matrix R[`r', 7] = r(table)[5,1] //R[`r', 4] - invnormal(0.975) * R[`r', 5] //CI lower
			matrix R[`r', 8] = r(table)[6,1] //R[`r', 4] + invnormal(0.975) * R[`r', 5] // CI upper
		 
			matrix list R
		 
			local r = `r' + 1 
		}

clear
svmat R


//local  cer=R[1,1]
//local  mse=R[2,1]

twoway (rcap R7 R8 R1, lcolor(navy)) ///
	(scatter R4 R1, mcolor(cranberry) yline(0, lcolor(black) lpattern(dash))), ///
	graphregion(color(white)) /// 
	xlabel( `cer'   "CER" `mse'  "MSE"  1.5 "1.5"  0 "0" `cer2' "2*CER"  `mse2' "2*MSE" 3 "3", angle(45))  /// 
	ylabel(-15(5)10, angle(horizontal)) legend(off) xtitle(" ") /*xtitle("Bandwidth") ytitle("RD estimate")*/  /// 
	title(" `label_`y'', All Schools, `label_`i'' ", bmargin(0)) //  

graph export "$output_dir/graph_bwsensitivity_8b4_`y'_post_`i'.png", as(png) replace


***

use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat "

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

}
}



/*
	Regression RD - PRIVATE
*/


use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat "

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

keep if sch_municip_pre==0
	
foreach i in 25 10{

local label_`i': variable label st_ind`i'

foreach y of global depvars{
	
		local label_`y': variable label `y'_post

		matrix define R = J(6, 9, .) // ( #of bw , #of statistics, missing )


		sum bw_cer_`y'_`i'_pv 
		local cer=round(r(min),0.1) // r(min) // 
		di `cer'

		sum bw_mse_`y'_`i'_pv
		local mse=round(r(min),0.1) // r(min) // round(r(min),0.1)
		di `mse'

		local cer2=`cer'*2
		local mse2=`mse'*2
		
		local r = 1
		foreach bw in `cer'  `mse' `cer2'  `mse2'   1.5  3 {
				
		
		reg `y'_post treat`i' st_ind`i' st_ind`i'_sq inter`i' inter`i'_sq  ${indepvars} `y'_pre  if abs(st_ind`i')<= `bw' , vce(cluster rbd)

			matrix list r(table)

			matrix R[`r', 1] =`bw' //round(`bw', 0.1) // name of bandwidth 
			matrix R[`r', 4] = r(table)[1,1] // rd estimate ****** used
			matrix R[`r', 5] = r(table)[2,1] //se_ of rd estimate ****** used
			matrix R[`r', 6] = r(table)[4,1] //2 * normal(-abs(R[`r', 4] / R[`r', 5])) // p-value
			matrix R[`r', 7] = r(table)[5,1] //R[`r', 4] - invnormal(0.975) * R[`r', 5] //CI lower
			matrix R[`r', 8] = r(table)[6,1] //R[`r', 4] + invnormal(0.975) * R[`r', 5] // CI upper
		 
			matrix list R
		 
			local r = `r' + 1 
		}

clear
svmat R


//local  cer=R[1,1]
//local  mse=R[2,1]

twoway (rcap R7 R8 R1, lcolor(navy)) ///
	(scatter R4 R1, mcolor(cranberry) yline(0, lcolor(black) lpattern(dash))), ///
	graphregion(color(white)) /// 
	xlabel( `cer'   "CER" `mse'  "MSE"  1.5 "1.5"  0 "0" `cer2' "2*CER"  `mse2' "2*MSE"  3 "3" , angle(45))  /// 
	ylabel(-15(5)10, angle(horizontal)) legend(off) xtitle(" ") /*xtitle("Bandwidth") ytitle("RD estimate")*/  /// 
	title(" `label_`y'', Private Schools, `label_`i'' ", bmargin(0)) //  
	
graph export "$output_dir/graph_bwsensitivity_8b4_`y'_post_`i'_pv.png", as(png) replace

***

use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat " // soc_post 

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

keep if sch_municip_pre==0

}
}



/*
	Regression RD - MUNICIPAL
*/


use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat "  

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

keep if sch_municip_pre==1 
	
foreach i in 25 10{

local label_`i': variable label st_ind`i'

foreach y of global depvars{
	
		local label_`y': variable label `y'_post

		matrix define R = J(6, 9, .) // ( #of bw , #of statistics, missing )


		sum bw_cer_`y'_`i'_mv 
		local cer=round(r(min),0.1) // r(min) // 
		di `cer'

		sum bw_mse_`y'_`i'_mv
		local mse=round(r(min),0.1) // r(min) // round(r(min),0.1)
		di `mse'


		local cer2=`cer'*2
		local mse2=`mse'*2
		
		local r = 1
		foreach bw in `cer'  `mse' `cer2'  `mse2'   1.5  3  {
				
		reg `y'_post  treat`i' st_ind`i' st_ind`i'_sq inter`i' inter`i'_sq  ${indepvars} `y'_pre if abs(st_ind`i')<= `bw' , vce(cluster rbd)

			matrix list r(table)

			matrix R[`r', 1] =`bw' //round(`bw', 0.1) // name of bandwidth 
			matrix R[`r', 4] = r(table)[1,1] // rd estimate ****** used
			matrix R[`r', 5] = r(table)[2,1] //se_ of rd estimate ****** used
			matrix R[`r', 6] = r(table)[4,1] //2 * normal(-abs(R[`r', 4] / R[`r', 5])) // p-value
			matrix R[`r', 7] = r(table)[5,1] //R[`r', 4] - invnormal(0.975) * R[`r', 5] //CI lower
			matrix R[`r', 8] = r(table)[6,1] //R[`r', 4] + invnormal(0.975) * R[`r', 5] // CI upper
		 
			matrix list R
		 
			local r = `r' + 1 
		}

clear
svmat R


//local  cer=R[1,1]
//local  mse=R[2,1]

twoway (rcap R7 R8 R1, lcolor(navy)) ///
	(scatter R4 R1, mcolor(cranberry) yline(0, lcolor(black) lpattern(dash))), ///
	graphregion(color(white)) /// 
	xlabel( `cer'   "CER" `mse'  "MSE"  1.5 "1.5"  0 "0" `cer2' "2*CER"  `mse2' "2*MSE"  3 "3" , angle(45))  /// 
	ylabel(-15(5)10, angle(horizontal)) legend(off) xtitle(" ") /*xtitle("Bandwidth") ytitle("RD estimate")*/  /// 
	title(" `label_`y'', Municipal Schools, `label_`i'' ", bmargin(0)) //  
	
graph export "$output_dir/graph_bwsensitivity_8b4_`y'_post_`i'_mv.png", as(png) replace


***

use "$data_dir/temp_b8_bw_nocovs.dta", clear

global indepvars "sch_municip_pre rural_rbd sch_stno_pre i.tipo_est_num i.mom_edu_pre i.dad_edu_pre i.pad_inc_pre i.sch_gse_pre i.hg i.sned_yr"  
global depvars "mate lect nat "  

la var mate_post "Grade 8 Mathematics"
la var lect_post "Grade 8 Reading"
la var nat_post  "Grade 8 Natural Sciences"

keep if sch_municip_pre==1 

}
}

erase "$data_dir/temp_b8_bw_nocovs.dta"





